Mesoscopic Model for Diffusion-Influenced Reaction Dynamics 
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A hybrid mesoscopic multi-particle collision model is used to study diffusion-influenced reac- 
tion kinetics. The mesoscopic particle dynamics conserves mass, momentum and energy so that 
hydrodynamic effects are fully taken into account. Reactive and non-reactive interactions with 
catalytic solute particles are described by full molecular dynamics. Results are presented for large- 
scale, three-dimensional simulations to study the influence of diffusion on the rate constants of the 
A + C ^ B + C reaction. In the limit of a dilute solution of catalytic C particles, the simulation 
results are compared with diffusion equation approaches for both the irreversible and reversible re- 
action cases. Simulation results for systems where the volume fraction <j> of catalytic spheres is high 
are also presented, and collective interactions among reactions on catalytic spheres that introduce 
volume fraction dependence in the rate constants are studied. 

PACS numbers: 02.70.Ns, 05.10.Gg, 05.20.Dd 



I. INTRODUCTION 

The dynamics of large complex systems often occurs 
on disparate time and space scales. The direct molecular 
dynamics simulation of the equations of motion for such 
systems is difficult because of this scale separation and 
the large numbers of molecules such systems may contain. 
Consequently, mesoscopic models play an important role 
in investigations of the dynamics of these systems. 

The use of Langevin and Fokker-Planck equations for 
Brownian motion is well known p], Q and these mod- 
els have been used in much wider contexts; for example, 
in investigations of reaction dynamics in the condensed 
phase yJl • Such stochastic models are useful when it is 
impossible or inappropriate to simulate the full dynamics 
of the system, including all solvent degrees of freedom. 

Suspensions of colloidal particles are also often treated 
using mesoscopic models of various types. While the dy- 
namics of the colloidal particles may be accurately mod- 
elled using Langevin dynamics, hydrodynamic interac- 
tions play an important role in dense colloidal suspen- 
sions. The friction tensors that enter the Langevin equa- 
tions depend on the colloidal particle configuration. To 
compute the frictional properties of dense suspensions, 
the intervening solvent is often approximated by the con- 
tinuum equations of hydrodynamics to determine the hy- 
drodynamic interactions among the colloidal particles. 

Other approaches for constructing mesoscopic dynam- 
ics of complex systems include the construction of ef- 
fective solvent models to be used in the context of full 
molecular dynamics simulations. 0] Such models allow 
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one to investigate systems of high complexity that can- 
not be studied by straightforward molecular dynamics 
simulation schemes. 

In this article show how diffusion-influenced reactions 
can be studied using a multi-particle mesoscopic dynam- 
ics. 0, Q In this dynamical scheme, particle positions 
and velocities are continuous variables and the dynamics 
consists of free streaming and multi-particle collisions. 
Multi-particle collisions are carried out by partitioning 
the system into cells and performing a specific type of 
random rotation of the particle velocities in each cell that 
conserves mass, momentum and energy. The hydrody- 
namic equations are obtained on long distance and time 
scales |6j and the model permits efficient simulation of 
hydrodynamic flows 0, Q • Since the dynamics is carried 
out at the particle level, it is straightforward to construct 
hybrid schemes where solute molecules that undergo full 
molecular dynamics are embedded in the mesoscopic sol- 
vent. 8] Hydrodynamic interactions among solute parti- 
cle are automatically accounted for in the multi-particle 
mesoscopic dynamics. ;9j The method has been gener- 
alized to treat phase segregating fluids with surfactants. 

m 

Diffusion-influenced reaction dynamics is widely used 
to model processes like enzymatic turnover or collision- 
induced isomerization in complex systems. Smolu- 
chowski constructed a continuum theory for such re- 
actions based on a solution of the diffusion equation. 
[Tl| In this article we focus on the reversible A + C ^ 
B + C reaction where a considerable body of research has 
concerned the development of refined theoretical mod- 
els FlI HI D3, 13 111 13 III Simulation schemes 
[3 l2Ct l2ll l22fl for three-dimensional diffusive reaction 
dynamics have been constructed. Diffusion-influenced 
reactions taking place in a dense field of catalytic parti- 
cles are strongly affected by perturbations of the diffusion 
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field arisin g fro m reactions at the different catalytic sites. 
IllHIlllllllllllllll This effect is similar to the hy- 
drodynamic interactions that enter colloidal suspension 
dynamics. We show how collective effects on diffusion- 
influenced reaction dynamics can be studied by simula- 
tions of a mesoscopic model for these systems. The meso- 
scopic multi-particle collision model allows us to simulate 
systems with tens of millions of particles for long times 
in order to determine power law decays and non-analytic 
catalytic particle density effects on the reaction rates. 

The outline of the paper is as follows. Section [n] 
sketches the mesoscopic multi-particle collision model 
and presents its generalization to multi-component sys- 
tems. The evolution equations that encode the multi- 
particle mesoscopic dynamics are presented in Sec. IIIII 
The computation of the diffusion coefficient, a necessary 
ingredient for the analysis of reaction dynamics, is given 
in Sec. IIVI In Sec. E] we show how the model can be 
generalized to treat chemical reactions. In particular we 
study the reaction A + C f± B + C that occurs upon 
collision with catalytic C particles. The simulation al- 
gorithms and simulation results for dilute and concen- 
trated suspensions of catalytic spheres are presented in 
Sec. I VII The conclusions of the investigation are con- 
tained in Sec. IVIII 



II. MULTI-COMPONENT MESOSCOPIC 
MULTI-PARTICLE DYNAMICS 

The mesoscopic dynamics we consider comprises two 
steps: multi-particle collisions among the particles and 
free streaming between collisions. 6] Suppose the system 
contains N particles with positions and velocities given 
by (X (Ar) ,VW) = (x 1 ,...,x 7V ,v 1 ,...,v w ). While the 
particle positions and velocities are continuous variables, 
for the purpose of effecting collisions, the system is di- 
vided into L cells labelled by the index £. Collisions occur 
locally in the cells in the following way: Rotation opera- 
tors Cj, chosen randomly from a set of rotation operators 
VL = {Cd\, . . . , <2ik} are assigned to each cell £ of the sys- 
tem. If a cell £ contains particles at time t and the cen- 
ter of mass velocity in the cell is Vj(t) = n^ 1 Y^iili v iMj 
the post-collision values of the velocities of the particles 
in the cell, v*, are computed by rotating the particle 
velocities relative to and adding to the result, 

v*(t)=V e (*)+^( Vi (t)-V e (t)) . (1) 

After the collision events in each cell, the particles free 
stream to their new positions at time t + r, 

x i (t + r)=x i (i)+vJ(t)T. (2) 

This simple dynamics has been shown to conserve mass, 
momentum and energy. The exact hydrodynamic equa- 
tions are obtained on macroscopic scales, and the sys- 
tem relaxes to an equilibrium Boltzmann distribution 



of velocities. |6( Consequently, the dynamics, although 
highly idealized, has correct behavior on macroscopic 
scales which are long compared to the effective collision 
times in the model. Since the dynamics is described at 
the particle level it is a simple matter to couple this 
mesoscopic dynamics to full molecular dynamics of so- 
lute species embedded in it. H |3(j The model is similar 
in spirit to Direct Simulation Monte Carlo |3l| but with a 
different discrete-time collision dynamics that simplifies 
the simulations and makes them more efficient. 

The mesoscopic dynamics for a multi-component sys- 
tem can be carried out in a similar way by generalizing 
the multi-particle collision rule. Suppose the TV-particle 
system comprises different species a = A, B, . . . with 
masses m a . In this case it is useful to introduce an opera- 
tor Qf that characterizes the species a of a given particle 
i. These operators have the following properties: 

e? ef = s aa , ; (3) 

i.e., particle i cannot be of different species at the same 
time; also, 

E 9 " = 1 ' ( 4 ) 

a 

so that particle i has to have some species type. The 
number of particles of species a is given by 

iV 

^-E "- ( 5 ) 

i=l 

There are many ways in which the multi-particle col- 
lision rule can be generalized for systems with several 
species and we consider one version that is consistent 
with the requirements that mass, momentum and energy 
be conserved. Let be the center of mass velocity of 
particles of species a that are in the cell £ at time i, 

v S a, (*) = -^7 £ e ^*)' (6) 

n 5 W i|xev 

where nl is the number of particles of the species a 
in cell £ with volume V at time t. The center of mass 
velocity of all n^(t) — n [ (*) particles in the cell £ at 

a 

time t is given by 

Vdt) = ^—1^ • (7) 

a 

In the model we adopt, two different types of multi- 
particle collisions occur. The first is a collision that in- 
volves particles of all species. To perform this collision, 
we use a rotation operator Cj which is applied to every 
particle in a cell as for single component system. The sec- 
ond type of multi-particle collision involves only particles 
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of the same species. The rotation operator uj a effects this 
collision and is applied to each particle of species a in the 
cell. Not only does it change from cell to cell and with 
time like Cj, but it also changes from species to species. 

The multi-particle collision process can be divided into 
these two independent steps. For the set of particles that 
are in the cell £, first we perform the all-species collision 
as 



V e + w c (vi-V, 



a > 



(8) 



where Vj the pre-collision velocity of the particle i and 
v'/ is the velocity after this step. Second, we apply the 
one-species rotation operator 



5>? (vf Q >+(2£(v?-Y 



11(a) 



(9) 



where V^°^ is the center of mass velocity of particles of 
species a after the all-species collision step. Note that 
u>£ is applied to all particles in the cell, but the u>^ are 
applied only on particles of species a. 

From Eqs. © and © the post-collision velocity of a 
particle may be expressed as 



v* = v £ +a, 5 (e?vf>-v 4 



e 



V 



(o) 



(10) 



III. EVOLUTION EQUATIONS 



The dynamics described above can be encoded in an 
evolution equation for the phase space probability den- 
sity, 

P( V W xW+vWr.i + r) 

= e^P(VWxW,t + T) 

= cp(y iN \x (N \t) , (ii) 

where the free streaming Liouville operator is, 

N 



(12) 



and N — N a is the total number of particles in the 

system. If we choose the rotation operators cD and Cj a 
randomly from the set f£, the collision operator may be 
written as, 

CP(vW,xW,t) = f i r ^|dV'W 

N 

a i—1 

-^(V^-V^-^K-V^)) , (13) 



where L is the number of cells. 

We may write the evolution equation in continuous 
time by introducing a delta function collision term which 
accounts for the fact that the multi-particle collisions oc- 
cur at discrete time intervals. We have 
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P(vW,xW,i) = (-L„ + C)P(VW,XW,() , 



.(14) 

where the collision operator C acts on the velocities of 
the particles at discrete times mr, and is defined as 

oo 
m=0 

(15) 

If Eq. I)14|l is integrated over a time interval mr—e to (m+ 
l)r— e we recover Eq. Ijll|l corresponding to multi-particle 
collision followed by free streaming. Instead, integration 
over the jump at t — (m+l)r yields an analogous discrete 
time equation with free streaming followed by collision. 

Assuming that the system is ergodic, then, in view 
of the conservation of mass, momentum and energy, the 
stationary distribution of the Markov chain in Eq. I|ll|) 
is given by the microcanonical ensemble expression, 

Po(v w xW) = W^f^e^Kii 2 -^ 

x*(EE 9 > Q (Vi-v)j , (16) 

\i=l a J 

where v is the mean velocity of the system, d is the di- 
mension and J\f is a normalization constant. If we inte- 
grate Pq over the phase space of all particles except par- 
ticle i, we obtain the Maxwell-Boltzmann distribution in 
the limit of large N. 

Figure ^ shows the results of a simulation of the ve- 
locity probability distribution for a system with volume 
V = 100 3 cells of unit length and N — 10 7 particles. 
The particles were initially uniformly distributed in the 
volume V and all particles had the same speed |v| = 1 
but different random directions. To obtain the results 
in this figure we assumed that the species were mechan- 
ically identical with mass m = 1 and used the multi- 
particle collision rule in Eq. l|lUfl with rotations ui£ and 
u>j? selected from the set il = {tt/2, —tt/2} about axes 
whose directions were chosen uniformly on the surface 
of a sphere (±7r/2 collision rule). This version of the 
collision rule for mechanically identical particles will be 
used in all calculations presented in this paper. The fig- 
ure compares the histogram of the x-component of the 
velocity with the Maxwell-Boltzmann distribution, 
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2tt 



(17) 



where (3 — {kgT)^ 1 , and confirms that this initial dis- 
tribution evolves to the Maxwell-Boltzmann distribution 
under under the mesoscopic dynamics. 
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Figure 1: Comparison of the simulated velocity distribution 
(histogram) with the Maxwell-Boltzmann distribution func- 
tion (solid line) for ksT =1/3. 



We may also write an evolution equation for any dy- 
namical variable a(V^ N ', X™) as, 

|a(vW XW i) = (L + C)a(vW,XW,t) , (18) 

where C has the same form as C in Eq. (|15(l with C re- 
placed by C, 



Ca(VW XW t) = ^£ jdV'W 

N 

xa(V'W,xW,t)niI ^(^- V « 

-^(V^-V^-^Vi-V^)) . (19) 

This equation is the starting point for the generalization 
to reacting systems in Sec. |V] 

IV. DIFFUSION 

A knowledge of the value of the diffusion coefficient 
is essential for the analysis of diffusion-influenced reac- 
tion kinetics. In this section we determine the diffusion 
coefficient as a function of the density from simulations 
of the mesoscopic multi-particle dynamics and derive an 
approximate analytical expression for its value. 

The diffusion coefficient is given by the time integral 
of the velocity correlation function. For the discrete time 
dynamics of the model, the time integral is replaced by 
its trapezoidal rule approximation, as shown by a discrete 
time Green-Kubo analysis, ja, H(j Thus, the diffusion 
coefficient D is given by 



D = ~{v x v x ) + ^2(v x v x (£t)) , 
1 t=i 



(20) 



where v x is the x-component of the velocity of a tagged 
particle in the system. (We suppress the species index 
a for the case of mechanically identical particles since 
all species have the same diffusion coefficient.) We have 
computed D using this expression as well as the formula 
for D in terms of the mean square displacement as a 
function of the mean particle density per cell p. The 
results are shown in Fig. 
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Figure 2: Comparison of the simulated diffusion coefficient 
(0) with the Boltzmann value (solid line). The ±7r/2 collision 
rule was used to obtain the results. The volume was V = 100 3 
and the temperature was ksT — 1/3. 

An approximate expression for D can be derived by 
assuming a single relaxation time approximation. If we 
suppose the decay is given by a single relaxation time, 
we have 



S ( PD ) 4 . (21) 



(v x v x (It)) ^ 
(v x v x ) V (v x v x ) 
The diffusion coefficient is then approximately given by 



D w --{v x v x ) + (v x v x ) ^ 



I'D 



2(1-7-23) 



The relaxation rate may be computed in the Boltz- 
mann approximation [fj], 



(vi x vix(t)) 



J ,,, ™— i ii ii J 



«Kv- Vl )I]>(vi)E'' 



jx i 



(23) 



3=1 



where, v\ x is the x component velocity of the single parti- 
cle 1. Since cross correlations between different particles 
are not present for self diffusion, we have, 



(vi x vi x (t)} 



1 1 UJ n=l ' i=l 



(24) 
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The x component of the post-collision velocity v* x may 
be written using Eq. (|10|l for the ±ir/2 collision rule dis- 
cussed above as 

»ix = ^ / ^( v * + [* ■ ( V(Q) - v )] ) 
dh {a) [ dh fii a) 



+ 



(4tt) 2 

x(n (Q) -n) n- (vi - V (q) ) 



(25) 



where, fi and are the normal vectors associated with 
the rotation operators uj and cD Q , respectively. As a result 
of this integration we obtain, 



(26) 



Assuming that particles of different species have the same 
mass, and substituting Eq. (|26|) into Eq. I)24[l. we find, 



rfVl^lj;) 



2 ( " + 2 

3n 



(vi) 
(27) 



For large enough p, we may approximate this expression 
by 



(wix^ix(r)) 
which yields, 



E 



2(l-e~P)+p 
3p 



(2 + n) 



(28) 



(29) 



Substituting ro in Eq. I)22|) the expression for the diffu- 
sion coefficient is 



D 



k b T (2p + l- e 



2m \ p — 1 + e 



(30) 



This analytic formula is compared with the simulation 
results in Fig. |3 where it is seen that it provides an ex- 
cellent approximation to the simulation results over all 
of the physically interesting density range. 



V. REACTIVE DYNAMICS 



The A and B particles undergo both non-reactive and re- 
active collisions with C, and the multi-particle collisions 
described in Sec. HTl among themselves. The macroscopic 
mass action rate law may be written as, 



df 



5nA(t) = — (kf + k r )SriA{t) = — kSnA{t) 



(32) 



where SfiA{t} = fiA(t)—n A q is the deviation of mean num- 
ber density of A particles from its equilibrium value, and 
k = kf + k r is the reciprocal of the chemical relaxation 
time. We have incorporated the fixed number density of 
the catalytic C particles into the rate constants. 

The microscopic evolution equation for this system 
may be written by simply augmenting the free streaming 
evolution operator in Eq. (|12|l with a Liouville operator 
L that describes the interactions of the A and B particles 
with the C particles. If the interactions of A and B with 
C are through continuous potentials, L takes the stan- 
dard form, L = F • Vp , where F is the force between the 
A and B particles and C and P is the vector of the the 
momenta of the particles. 

For the purposes of calculation and illustration, we 
adopt a model where the C particles are fixed in space 
and have radius a. The A and B particles either bounce 
back from the catalytic spheres without changing their 
identity or react with probability pr. In this case the 
evolution equation for any dynamical variable in the sys- 
tem is given by 

|a(xW VW *) = (L ± L ± + C)a(xW V^.t) , 
at 

(33) 

where the ± signs apply for t > and t < 0, respec- 
tively. The Liouville operators L± describing the reactive 
and non-reactive collisions with the catalytic particles are 
given by 

M N 

a j — 1 i—1 

M N 

a j—1 i—1 



xbi 



lie? 



(34) 



Here f y = (xj — xy)/ is a unit vector along the line 
of centers between particle i and the catalytic sphere j, 
Tij = |x, — Xj| is the magnitude of this vector and the 
operator 6^ converts the velocity of particle i to its post- 
collision value after collision with the catalytic sphere j, 



Next, we consider a reactive system with M finite-sized 
catalytic spherical particles (C), and a total of N — Na + 
Nb A and B particles which react with the C particles 
through the reactions, 

A + C%B + C . (31) 



&*i(vi,V 2 , ...,Vi,.. . ,Vjv) = (vi,V 2 , . . . ,V*, . . . ,Vjv) ■ 

(35) 

For bounce-back dynamics we have v* = — Vj. The op- 
erator V aa acts on the species labels to effect reactive 
collisions so that V aa '&? = 8f' where a' = B if a = A 
and vice versa. The factor 7 accounts for the possibility 
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that the forward and reverse reactions occur with dif- 
ferent probabilities leading to an equilibrium constant 
K eq = 7 _1 which is different from unity. 



VI. RESULTS 
A. Simulation method 



Rate law 

The chemical rate law for this system my be derived 
by taking the dynamical variable a to be the deviation 
of the number of particles of species A from its average 
value, x — N A — < N A >= SNa — —SNb, where 



N 



n a = J2®1 



(36) 



The angular brackets < • • • > signify an average over 
an equilibrium ensemble where the numbers of A and B 
molecules fluctuate but their sum is fixed, Na + Nb = N 
Starting with Eq. for t > and using standard pro- 
jection operator methods |3^ | we may write a generalized 
Langevin equation for \(t) in the form, 



d u\ f m < L - X)X > M 

-7lX{t) = fxW —X(t) 

dt <XX> 



f 

Jo 



dt x{t 

<XX> 



f ) , (37) 



where we have introduced the projection operator Va =< 
ax >< XX >~ 1 X an d its complement Q = 1 — V. The 
random force is f x (t) = exp[QL + t]QL + x- 

Averaging this equation over a non-equilibrium ensem- 
ble where x does not fluctuate yields the generalized 
chemical rate law, 

—Sn A (t) = —5n A (t) 

dt <XX> 



f 

Jo 



^<(L-x)e«M' QL+x> 

< XX > 



t'). 



The contribution 

< (L-X)X > 
< XX > 



= k of (l + K- 1 ) , 



(38) 



(39) 



determines the initial rate arising from direct collisions 
of the A and B particles with the catalytic spheres. For 
bounce-back collision dynamics of the A and B species 
with the catalytic sphere C, we have 



2 /87rfc B T\i/2 

«o/ = Pro- nc 

\ m J 



(40) 



where nc is the constant number density of catalytic 
spheres. The memory term accounts for all diffusion- 
influenced effects arising from recollisions with the cat- 
alytic spheres. 



The simulation of model is carried out in a cubic box 
with sides Lb and periodic boundary conditions. The 
centers of the spheres of radius a are located in this box, 
taking care to preserve periodic conditions on the edges 
when the spheres lie partially outside the cube. Once the 
catalytic spheres are placed in the box, the initial posi- 
tions of the particles are assigned values that are within 
the cube but outside the spheres. The velocities are cho- 
sen from a Maxwell-Boltzmann distribution. 

Given the initial distribution of particles and particle 
velocities, the simulation begins by grouping the parti- 
cles in cubic cells of size 1 within which the multi-particle 
collision operators act to change the velocities of all par- 
ticles, preserving their positions. Then the displacement 
of each particle is computed using the post-collision ve- 
locity, taking into account the periodic boundary condi- 
tions of the cube and the bound-back collisions with the 
spheres. When a particle hits a sphere it may react with 
probability pa, and the sign of its velocity is changed. 
Collisions between particles and spheres occur in contin- 
uous time in the interval [t, t + r]. When many catalytic 
spheres are present a particle may hit several spheres in 
one unit time r. 

Once all the particles have been moved, the time ad- 
vances one unit r and the particles are regrouped to apply 
the multi-particle collision rule again. 



B. Single catalytic sphere 

In order to test the utility of the mesoscopic model we 
investigate a system that contains a dilute distribution of 
independent catalytic C particles so that the dynamics 
may be described by considering a single C particle (la- 
belled 1) with radius a in a medium of A and B particles. 
In the case where A particles are converted irreversibly 
to B upon collision with C the chemical rate law takes 
the form, dn A if)/dt — —kf(t)nA(t), where is the 

time dependent rate coefficient. If the dynamics of the 
A density field may be described by a diffusion equation, 
we have the standard partially absorbing sink problem 
first considered by Smoluchowski. To determine the 
rate constant we must solve the diffusion equation 



dn A (r,t) 
dt 



D A n A {r,t) , 



subject to the boundary condition |33l | 

4nDa 2 r ■ (VnA)(ro-, t) = kofn A (ra,t) 



(41) 



(42) 



This equation assumes that the continuum diffusion 
equation is valid up to a > a, which accounts for the pres- 
ence of a boundary layer in the vicinity of the the sphere 
surface where the continuum diffusion description should 
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fail. The resulting expression for the time-dependent rate 
coefficient is Hj 



150 



*/(*) = 



c / + ki 

k 2 
K o/ 



K f + 

x crfc 



ko£\fDt 
fcn/U 2 



1 



(43) 



Here fcc = 47r<j.D is the rate constant for a diffusion 
controlled reaction for a perfectly absorbing sphere. 

The time-dependent rate coefficient kf(t) may be 
determined directly from the simulation by moni- 
toring the A species density field and computing 
~(dnA(t)/dt)/nA(t). The results of such a computa- 
tion for irreversible reaction (7 = 0) with probability 
Pr = 0.5 is shown in Fig. [3] The system size is 100 3 




1000 



Figure 3: Plot of the time dependent rate constant kf(t)/nc 
versus t for a = 10. Fhe solid line is theoretical value of kf(t) 
using Eq. 1431 and a = a + 1. 



kf 75 - 




a 



300 



kf 150 - 




a 



Figure 4: Plot of kf/nc (©) versus a, the radius of the 
catalytic sphere. The initial value kf(t = 0) = kof (□) 
is also plotted versus a in this figure. The solid lines are 
the theoretical values of these quantities determined from 
k^ 1 — (fco/(l + Keq))^ 1 + fc^, 1 . (a) Irreversible reaction 



{K~ = 0) with PR 
with pr — 1. 



0.5. (b) Reversible reaction (K' e 



volume units and there is a sphere of radius a = 10 lo- 
cated in the center of the system. The simulation starts 
with N(0) = N A (0) = 10 7 particles of species A with 
unit mass uniformly distributed in the space. The initial 
velocities are Maxwell distributed with /esT/m = 1/3. 
The time dependent rate coefficient starts at kof and de- 
cays to its asymptotic value kf. In our mesoscopic model 
the continuum theory cannot apply on the scale of one 
multi-particle collision cell, so we have taken a = a + 1 
to approximately account for the microscopic boundary 
layer. One sees good agreement between the simulation 
and diffusion theory results. 

In Fig.0Ji we plot the values of kf extracted from the 
simulation data in this way versus the radius of the cat- 
alytic sphere. The figure shows the increasing importance 
of diffusion-influenced effects on the value of the rate con- 
stant as a increases. While fco/ grows quadratically with 
a in accord with Eq. 140fl . we see that kf grows more 



slowly and approaches the diffusion- limited value of ko, 
which depends linearly on a for large a. The theoretical 
estimate, kj 1 — k^ + fc^ 1 , is in good agreement with 
the simulation results. 

A similar calculation can be carried out for the re- 
versible case (7 = 1 andpji — 1). For reversible reactions 
the chemical relaxation rate k(t) is given by Eq. I|43l) with 
kof replaced by k$ = fco/+^Or = kof(l+K~ ) and, there- 
fore, fc _1 = k^ 1 + kp 1 . ^3 For our simulation conditions 
K eq = 1 so that ko — 2fco/- Also kf — k r . In Fig.0J) we 
plot the simulation values of kf for the reversible reaction 
and compare them with the diffusion equation formula. 
Once again good agreement is found. The effects of dif- 
fusion appear at somewhat smaller values of a since ko is 
larger for the reversible reaction and the diffusion-limited 
value of the rate constant is reached at smaller values of 
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C. Random distribution of catalytic spheres 

If instead of a single catalytic sphere we have a ran- 
dom distribution of M spheres of radius a in the vol- 
ume V, the rate constant will depend in a non-trivial 
way on the catalytic sphere density or volume fraction 
4> = 4tt<t 3 M I (3V) . The reactions at one sphere surface 
will alter the A and B particle density fields there. From 
the perspective of a continuum diffusion equation ap- 
proach, since the diffusion Green function which couples 
the dynamics at the different spheres is long ranged, the 
interactions from many catalytic spheres determine the 
value of the rate constant. The problem is analogous to 
the long range interactions that determine hydrodynamic 
effects on the many-particle friction coefficient. There 
have been a number of studies of the volume fraction de- 
pendence of the rate constant [iH El El El El El |28| . 
These derivations rely on resummations of classes of 
interactions among the reacting spheres or other tech- 
niques. 

The chemical relaxation rate for a system with a ran- 
dom distribution of catal ytic spheres with volume frac- 
tion 4> is given by [H Ell2Hl 



k(4>) = k 



(k f + fc 0r ) £ 
(kof + k 0r 



r30 



1/2 



(44) 



k D x . The first 



where, as earlier, k = (fco/ + ko r )~ 
finite density correction to the rate constant depends on 
the square root of the volume fraction. This non-analytic 
volume fraction dependence arises from the fact that the 
diffusion Green function acts like a screened Coulomb po- 
tential coupling the diffusion fields around the catalytic 
spheres. As in the Debye theory of electrolytes, one must 
sum an infinite series of divergent terms to obtain the 
non-analytic <p dependence. 

The mesoscopic multi-particle collision dynamics fol- 
lows the motions of all of the reacting species and their 
interactions with the catalytic spheres. Consequently, all 
many-sphere collective effects are automatically incorpo- 
rated in the dynamics. We have carried out simulations 
of the chemical relaxation rate constant k(<fi) as a func- 
tion of the volume fraction of the catalytic spheres for a 
reversible reaction with 7=1 {K~g = 1) and p^ = 0.25 
as well as an irreversible reaction with 7 = [K~^ = 0) 
and pa = 0.5. For this choice of parameters the theoret- 
ical formula predicts that k((f>) for the reversible reaction 
is equal to kf((f>) for the irreversible reaction. Our simula- 
tions were performed for systems with a volume fraction 
<j) of catalytic spheres with radius a = 3 in a system of size 
100 3 multi-particle cells and an initial number density of 
A particles, n^(0) = 10 per cell. The results shown in 
Fig. [S] were obtained from an average over five realiza- 
tions of the random distribution of catalytic spheres. We 
see that the simulation results confirm the existence of a 




Figure 5: Relaxation rate coefficient k(<f>)/nc as a function 
of the square root of the volume fraction cjf> 1//2 for a = 3 and 
ksT = 1/3. Irreversible reaction k[{(j>) (•). Reversible reac- 
tion k(<f>) (0). For some values of (j> the two cases cannot be 
distinguished in the figure because the data points overlap. 
The solid line is determined using Eq. llHt . 



dependence on the volume fraction for small volume 
fractions. As predicted by the theory for the chosen pa- 
rameter values the reversible and irreversible data over- 
lap, even in the high volume fraction regime. For larger 
volume fractions the results deviate from the predictions 
of Eq. Ij44(l and the rate constant depends much more 
strongly on the volume fraction. In this regime the dif- 
fusion coefficient is also modified as a result of collisions 
with the catalytic spheres and this effect also contributes 
to the deviation. 

From these results we conclude that the mesoscopic 
multi-particle collision dynamics provides a powerful tool 
for the exploration of concentration effects on diffusion- 
influenced reaction kinetics. Such concentration depen- 
dence is often difficult to explore by other means. 



VII. CONCLUSION 

We have demonstrated that large-scale simulations of 
diffusion-influenced reaction kinetics are possible by us- 
ing the mesoscopic multi-particle collision model. With 
this model the dynamics of tens of millions of particles 
interacting with hundreds of catalytic spheres could be 
followed for long times to obtain the rate constants char- 
acterizing the population decay. Such simulations would 
be very costly using full molecular dynamics methods. 

Since the dynamics is followed at the (mesoscopic) par- 
ticle level, a number of noteworthy features of the dy- 
namical scheme are worth mentioning. From a technical 
point of view the dynamics is stable and no difficulties 
like those associated with discretizations of the diffusion 
equation or boundary conditions arise. Reversible and 



9 



irreversible reaction kinetics may be treated in similar 
fashion. All details of interactions arising from compe- 
tition among the catalytic spheres in a dense suspension 
are automatically taken into account; thus, screening ef- 
fects enter naturally in the dynamics. 

The model may be generalized to any reaction scheme 
and is not restricted to the simple A + C ^ B + C re- 
action with catalytic C particles discussed in this paper. 
Since solute molecules embedded in the mesoscopic sol- 
vent evolve by full molecular dynamics (without solvent- 
solvent interactions), the model will be most efficient 



when solvent-solvent interactions are a major time lim- 
iting factor in the simulation. This could be case for 
conformational changes of large molecules in solution, re- 
actions involving energy transfer in solution, etc. Thus, 
the model should find applicability in a variety of cir- 
cumstances when diffusion-influenced reaction kinetics is 
important. 
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